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We investigate the innermost stable circular orbit (ISCO) of a test particle moving on the equa- 
torial plane around rotating relativistic stars such as neutron stars. First, we derive approximate 
analytic formulas for the angular velocity and circumferential radius at the ISCO making use of an 
approximate relativistic solution which is characterized by arbitrary mass, spin, mass quadrupole, 
current octapole and mass 2 4 -pole moments. Then, we show that the analytic formulas are accurate 
enough by comparing them with numerical results, which are obtained by analyzing the vacuum 
exterior around numerically computed geometries for rotating stars of polytropic equation of state. 
We demonstrate that contribution of mass quadrupole moment for determining the angular veloc- 
ity and, in particular, the circumferential radius at the ISCO around a rapidly rotating star is as 
important as that of spin. 
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I. INTRODUCTION 

Observations of the low-mass X-ray binaries (LMXBs) 
with the Rossi X-ray Timing Explorer have revealed 
quasi-periodic oscillations (QPOs) around a frequency of 
~ 1 kHz 0. At present, more than 10 sources of kHz 
QPOs have been found One of the most impressive 
features of kHz QPOs is their very high frequency itself. 
The LMXBs are considered to be systems which include 
a neutron star of mass M ~ 1.4M Q , where Mq denotes 
the solar mass, and the accretion disk around the neutron 
star. The Kepler frequency of the accretion disk is 
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where R denotes the circumferential radius around the 
neutron star. If the origin of the kHz QPOs is certain 
oscillation frequencies of the accretion disk surrounding 
a neutron star of low magnetic field, they must be gener- 
ated at less than ten Schwarzschild radii of the neutron 
star. This means that the kHz QPOs may bring us a 
chance to explore general relativistic effects [|| . 

Several authors have recently suggested that at least 
some of the kHz QPOs may be related to the Kepler fre- 
quency at the innermost stable circular orbit (ISCO) of 
the accretion disk around a neutron star Jl]] . One of the 
most strong reasons is that in many sources, the maxi- 
mum frequency of the kHz QPOs is in narrow range be- 
tween 1.1 and 1.2kHz, although they are thought to have 
very different mass accretion rates and magnetic fields. 



Due to the fact that the ISCO is determined by the prop- 
erty of the central neutron star, but not by the properties 
of the accretion disk, such as the mass accretion rate, it 
has been suggested that the origin of the kHz QPOs of 
the highest frequencies may be the Kepler frequency at 
the ISCO. If this is true, it means that we have a great 
opportunity to investigate general relativistic effects 

I- 

Another remarkable feature of kHz QPOs is that they 
exhibit some evidences that in the center of their sources, 
rapidly rotating neutron stars are involved as follows: (1) 
many sources display two peaks of the kHz QPO, and 
the frequency difference between the twin peaks (A/) 
changes only slightly with time although the frequency 
of each peak changes jlj]; (2) some sources which pos- 
sess twin peaks also exhibit very coherent oscillations 
of several hundred Hz in X-ray bursts Jl[, and the fre- 
quencies change little with time during the bursts. Fur- 
thermore, they are approximately equal to or twice A/. 
Since the spin frequency of a neutron star is the only can- 
didate which changes only slightly on short time scales, 
the origin of the frequency difference between the twin 
peaks of QPOs and the oscillation frequencies in the X- 
ray bursts seem to correspond to the spin frequencies (or 
twice them) of neutron stars [jlj. This means that QPO 
sources include rapidly rotating neutron stars. 

Since the ISCO is determined by the geometry around 
the star, it is important to ask if the geometry around a 
neutron star can be approximated by that around a black 
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hole. If the electric charge is neutral, stationary black 
holes must be of the Kerr type due to the uniqueness the- 
orem Kerr black holes have mass, spin, quadrupole 
moment and so on, but multipole moments higher than 
the quadrupole are expressed in terms of the mass M and 
the spin angular momentum J = Si = qM 2 (\q\ < 1) 
as M 2l = M{iqM) 21 and S 2 i-i = -iMiiqM) 21 - 1 (I = 
1, 2, 3, -)P| |§, where M; and Si denote the mass and 
current moments, respectively. This means that the ge- 
ometry outside the black hole horizon is expressed only 
in terms of M and q. As a result, the ISCO is deter- 
mined solely by them. However, this is not the case for 
neutron stars. In the neutron star case, multipole mo- 
ments higher than and including the quadrupole do not 
depend on the mass and spin in such a simple manner, 
and they are determined when the equation of state of 
the neutron star is given. 

The purpose of this paper is to point out the signifi- 
cance of the multipole moments (in particular, the mass 
quadrupole moment) of rapidly rotating neutron stars 
in determining the ISCO around them. This is due to 
the following fact: In the case of a Kerr black hole, 
the magnitude of the quadrupole moment is denoted by 
M 2 = -q 2 M 3 [§, and \M 2 \ is always smaller than M 3 . 
In this case, the effect of the quadrupole moment is not 
important except in the case q ~ 1. However, in the 
case of a rotating neutron star, \M 2 /q 2 M 3 \ may be much 
larger than unity ~ 10, and hence for rapidly rotating 
neutron stars which seem to be located at the center 
of QPO sources, \M 2 \ may be larger than qM 3 for the 
case q > 0.1 ||]. In such a case, the effect due to the 
quadrupole moment is as important as that due to the 
spin. 

The paper is organized as follows. In Sec. II, we derive 
approximate analytic formulas for the angular velocity 
and circumferential radius of the ISCO around a rotat- 
ing object characterized by its mass, spin angular mo- 
mentum, mass quadrupole, current octapole, and mass 
2 4 -pole moments. In Sec. Ill, we perform numerical com- 
putations for obtaining stationary, axisymmetric space- 
times of rotating stars, and making use of the numerical 
results, we demonstrate that the accuracy of our formu- 



* In this paper, we use units of G = 1 = c, where G and c 
are the gravitational constant and speed of light. 



las derived in Sec. II are accurate. Sec. IV is devoted to 
a summary. 

II. APPROXIMATE ANALYTIC FORMULAS 
FOR THE ANGULAR VELOCITY AND 
CIRCUMFERENTIAL RADIUS AT ISCO 

A. Basic equations 

The line element of the vacuum exterior outside a sta- 
tionary, axisymmetric rotating object is generally written 



ds 2 = -F(dt - LodLpf 



1 r 

F 



e 2l {dp 2 +dz 2 )+ p 2 dip 2 . 

(2.1) 



Throughout this paper, we assume that the spacetime 
has reflection symmetry with respect to the equatorial 
plane z = 0, so that F, uj, and 7 are functions of p and 
\z\. 

Our purpose is to derive approximate formulas for the 
angular velocity and circumferential radius at the ISCO 
of a test particle on the equatorial plane around a rotat- 
ing object. In the case when the test particle stays on 
the equatorial plane (z = 0), geodesic equations can be 
integrated to give 



dt _ Eg w + £g ttp 
dr g 2 

dip _ Eg tlp + lg tt 
dr g 2 

l//ip 2 = _i + e 2 i££ + 2El 9 -^ +£ 29 -^ 



dr 



92 



92 



92 



-V(p), 



(2.2) 
(2.3) 

(2.4) 



where E and £ denote the specific energy and specific 
angular momentum of the test particle. g tt — —F, g ttp = 
Fu, g vv = -Fu 2 + p 2 /F, and g 2 = g 2 v - gttg vv (52 = p 2 
in the present line element). For circular orbits, relations 
of the angular velocity fi, E, and I are derived from the 
conditions V = = dV/dp as 



CI = 



d:f 



dt 

E = - 



9w,p 



gu + gtipCi 



£ = 



\J-gu - 2gt v Ci - g vv Ci 2 ' 

9 t v + ffyy^ 

\/-gtt - 2g tv fl - g vv Cl 2 



(2.5) 
(2.6) 
(2.7) 
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For simplicity, we only consider prograde orbits in this 
paper. A circular orbit is stable (unstable) if 



d 2 V 



1 



dp 2 g 2 \ dp 



d 2 92 _ E 2 d 2 g w 



dp 2 
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d 2 gt v 
dp 2 



n2 d2 9tt 



dp 2 



(2.8) 



is positive (negative) . Hence, the coordinate radius p and 
Q at the ISCO are determined from the condition where 
d 2 V/dp 2 is vanishing. 

Note that Eq. ( |2.8| ) is independent of metric function 
7. Also, the angular velocity £1 and the circumferential 
radius {y/g vv = R) are independent of 7. Thus, we only 
need F and u> in the following. 

B. Metric from the Ernst potential 

A stationary axisymmetric vacuum geometry is com- 
pletely determined by the Ernst potential ||, which is 
defined as 



£ = F + i^j 



^/p 2 + z 2 +i 



(2.9) 



where F — —gtt and £ is a complex potential. If we know 
ij), u> is calculated as 



r p' dip , 



constant z 



(2.10) 



Thus, once we know £, we have all the necessary infor- 
mation. 

£ has the property that it can be expanded as M Jul 



a 2j,k 



p 2 ^Z k 



j=0 k=0 



[p 2 + Z 2 ) 2 3+ k ' 



(2-11) 



where <2j,fc are complex constants in which information of 
the multipole moments of spacetime is completely con- 
tained. Note that a,j^ is non-zero only for non-negative, 
even j and non-negative k. Also, because of reflection 
symmetry with respect to the equatorial plane, a^fe is 
real for even k, and pure imaginary for odd k. Note that 
for investigation of the ISCO on the equatorial plane, we 
only need ctjfi and 

Fodor, Hoenselaers, and Perjes (FHP) Q show that 
all the components of a^fc are derived by the following 
recursive relation 



1 



&r,s+2 — 



(* + !)(* + 2) 



-(r + 2)V+2 



+ E a W a *-k-p,s-l-j x 
k,l,p,j 

Kj {p 2 + j 2 -Ap- 5j - 2pk - 2jl - 2) 
+a p+2 ,j-2(p + 2){p + 2-2k) 



+a P _2j+2(7+2)0' + l-20} 



(2-12) 



where the sum is taken for < < r, < / < s + 1, 
< p < r — fc, and — 1 < j < s — I |hJ , and a* k denotes 
the complex conjugate of As pointed out in ref. p0[ , 
it can be shown that a r ^ s +2 (s > 0) is a function of Oj.o 
and aj-i^i with j < r + s + 2. This means that if we 
know Oj.o and aj.i for j > 0, the entire spacetime metric 
(and of course, the metric on the equatorial plane) are 
completely determined. Note that if we know Oj.o up to 
j = 2n and a^i up to j = 2n — 2, we can calculate ao,fe up 
to k — 2n. Conversely, if we know ao,fc up to k — 2n, we 
can calculate djfi up to j = 2n and aj t \ up to j = 2n — 2. 
In principle, we can calculate terms of a /- up to arbi- 



trary large j and k using Eq. (2.12). In practice, however, 
we have to truncate higher terms. To access an appro- 
priate truncation point, we can make use of the solution 
for slowly rotating black holes of the mass M and the 
angular momentum J(= Si) -C M 2 . In this case, 



Mp 



izpj 



vV + M 2 ( P 2 + M 2 )V- 



+ 0{z 2 ), (2.13) 



where we have neglected terms of 0(J 2 /M 4 ). Note that 
p is related to the Schwarzschild radial coordinate r s as 
p = ^r s (r s — 2M). The Schwarzschild coordinate radius 
of the ISCO is r s = 6M(1 - ^ 8/27. J/M 2 ), so that the 
radius of the ISCO in p is given by 



Pisco =V24M-^. 



(2.14) 



We have investigated how the estimated values of pisco 
differ from the above one by expanding £ in terms of 
0(M/p) as 



M 2 
2p 
zJ 



3M 4 5M 6 35M 8 



63M 



10 



8p 4 16p 6 128p 8 256p 10 
3M 2 15Af 4 35Af 6 315M 8 



2p 2 



&P 4 



16p e 



128p 



-0{p 



(2.15) 



We have found that in the case when we take up to 
0(p~ 6 ), 0(p~ 8 ), and 0{p~ w ) terms, the errors in pisco 
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are about 4x 10 4 , 10 5 and 4x 10 7 , respectively, for the 
coefficient of M, and about 2 x 10~ 3 , 10~ 4 and 4 x 10~ 6 , 
respectively, for the coefficient of J/M. The order of 
magnitude of the errors in fiisco have turned out to be 
almost the same as those in pisco- Thus, we expect that 
this method will generate a fairly accurate approximate 
formula even if we include the higher multipole moments. 
In the following, we take a^o up to j — 10 and a^i up 
to j = 8, i.e., we take into account all the terms up to 
0(p~ w ) consistently. In other word, we calculate ao,fc up 
to k = 10, and neglect ao,fc of k > 11. 

C. Results 

Our strategy for determining the ISCO is as follows. 
First, we assume that the spacetime is characterized by 
mass M, spin angular momentum J(= S\) = qM 2 , mass 
quadrupole Af 2 — — Q 2 M 3 , current octapole moments 
S3 = -Q3-M 4 , and mass 2 4 -pole M 4 = Q±M b , and neglect 
the higher multipole moments. Note that q, Q 2 , 93, and 
Q4 are positive for a rotating star in a prograde spin. 

In the case of a rotating neutron star, we may assume 
q ~ O(ei), Q 2 ~ 0(e 2 ), 93 ~ 0(eie 2 ), Q± ~ 0(e|)> where 
£i, e 2 <C T because q and Q 2 are expected to be less than 
unity. For a slowly rotating neutron star, Q 2 ~ 0(q 2 ), 
so that e 2 ~ 0{e\). In this case, Q 2 ~ 0(e 2 ), 53 ~ 0(ef ), 
Q4 ~ O(ef). However, for a rapidly rotating neutron 
star, Q2 can be as large as q, so that whenever the terms 
proportional to q 2 make a significant contribution, we 
should also take into account the terms proportional to 
Q2, Qa, qq3, and so on. This is the reason why we take 
into account S3 and M4. 

Hereafter, we expand all the quantities by means of 
e = e\ by formally setting q — ► eq, Q 2 — > e 2 Q 2 ,<?3 — > 
e 3 (/3, and Q4 — > e 4 Q4, and retain all the terms up to 
0(e 4 ) consistently. Thus, the formulas derived below are 
accurate up to 0(e 4 ) for a slowly rotating neutron star. 
Even for a rapidly rotating neutron star, the formulas 
include all the terms of O(ei), 0(e 2 ), 0{t\), 0(e%), and 
0(eie 2 ). Hence, they are still accurate to 0{Q\). 

As mentioned in the previous subsection, we need ao,fc 
up to k = 10. Relation between multipole moments and 
ao,k for < k < 10 have been already given by FHP ||. 
Hence, by using them, we can calculate a,2jfl (1 < j < 5) 
and a 2 j,i (1 < j < 4) up to 0(e 4 ) using Eq. gl). 
(The explicit forms of ao.fe) a 2j,o and i2j,i are shown in 



Appendix A). Once £ is determined up to 0(p 
and u> are straightforwardly obtained from Eqs 



and (2.10) in the following form; 



F = 1 



11 



3=1 v 1 j 



11 



i=3 



+ o(p- 12 ), 



10 ), F 

dH) 

(2.16) 
(2.17) 



where Cpj and C WJ are functions of q, Q 2 , 53, and 
Q4. Substituting these equations into d 2 V/dp 2 — 0, we 



rewrite Eq. (2.8) in the form 



4 

= o, 

i=0 



(2.18) 



where the coefficients A; are independent of e but depend 
on q, Q2, (73, Q4, and p. These coefficients can be explic- 
itly written out by substituting Eqs. (2.16) and (2.17) 



into Eq. (making use of Eqs. (gj) and (gj)) and 

gathering in powers of e. We then look for a solution in 
the form 



i=0 



(2.19) 



where we know that c = V24M and c\ = — 10gM/3 



Imposing that Eq. ( 2.18 ) holds in each order of e, it is 
rewritten into three algebraic equations for c 2 , C3, and C4 
as 



1 d 2 A 

2 dp 2 ' 



dAn 



Id 3 A 



6 dp : 



dp 
dp 1 



c 2 



dAi 
dp 



Cl 



A 2 = 0, 



-ClC 2 



dA 
dp 



c 3 



1 d 2 Ai 2 dAi 
+ 2^ Cl + ^ C2 



■ci 



1 d 4 A 



24 dp 4 



dA 



ld 3 A 2 
2^ 3 - ClC2 + 
ld 3 Ai , 



dp 4 6 dp 3 



dA 2 
dp 

1 rf 2 A 2 

2 dp 2 

d 2 A! 



ld 2 A 2 



dA 2 



-c 2 



dp 2 
dA 3 



-ClC 2 



"Cl 



A3 = 0, 

2cic 3 ) 

dAi 
dp 

A 4 = 



(2.20) 



(2.21) 



C3 



(2.22) 



2 dp 2 1 dp ' dp 

where A^ and their derivatives are evaluated at p = cq. 
From these algebraic equations we obtain 

c 2 ~ (-1.41671<? 2 + 1.07648Q 2 )M, 

c 3 ~ (-1.45081q 3 + 1.60313£?Q 2 - 0.32561q 3 )M, 

c 4 ~ (-1.87941q 4 + 2.74953g 2 Q 2 - 0.37568Q 2 . 

+ 0.09377Q 4 - 0.70086 W3 )Af. (2.23) 



4 



Using these approximate solutions, we reach approximate 
expressions for the angular velocity and circumferential 
radius of the ISCO, ^isco and -Risco, as follows. 

Oisco = 6 Jq M (* + 0.74846 9 + 0.78059q 2 - 0.23429Q 2 

+ 0.98094g 3 - 0.64406gQ 2 + 0.07432 ?3 
+ 1.38118 9 4 - 1.41729g 2 Q 2 + 0.12798Q^ 
- 0.02129Q 4 + 0.25028^3) , (2.24) 

R ISCO = 6M(l ~ 0.54433<7 - 0.22619g 2 + 0.17989Q 2 

- 0.23002 9 3 + 0.26296gQ 2 - 0.05317g 3 

- 0.29693? 4 + 0.44546g 2 Q 2 - 0.06249Q2 

+ 0.01544Q 4 - 0.11310^3) • (2-25) 

For the case of Kerr metric, Q 2 = q 2 , q% = q 3 , and Q4 = 
g 4 , we obtain 

^isco = 1 f 1 + 0.74846? + 0.54630 9 2 + 0.41120g 3 
6v 6M v 

+ 0.32085g 4 ), (2.26) 
-Risco = 6M (l - 0.54433g- 0.04630g 2 - 0.02023g 3 

-0.01162g 4 ). (2.27) 
On the other hand, the exact Kerr solution gives fl2]j 

n*°™ = (l + 0.74846? + 0.54630g 2 + 0.41108g 3 

6v 6M \ 

+ 0.31991g 4 + O(< ? 5 )), (2.28) 
R fsco = 6M ( 1 - 0.54433g- 0.04630 g 2 - 0.02016g 3 

- 0.01H0? 4 + 0(g 5 )Y (2.29) 



Thus, the error of the coefficients in our approximate 
formulas of f^isco and -Risco is less than 10~~ 4 for 0(e 2 ) 
terms, less than 10~ 3 for 0(e 3 ) terms, and ~ 10 -2 for 
0(e 4 ) terms. Hence, for slowly rotating neutron stars 
(q < 0.1), our formulas are accurate enough. Even for 
very rapidly rotating neutron stars of q, Q 2 ~ 0.5, we 
may expect that they will yield very accurate values. 



We note that in Eqs. ( 2.24Q and (2.25), the signs of 
coefficients of the terms including Q 2 such as Q2, C1Q2 
and q 2 Qi are opposite to those of q k (k = 1 — 4) terms. 
As shown in ref. || , Q 2 /g 2 is always larger than unity for 
a rotating neutron star, and hence Q 2 may be of 0(q) for 
rapidly rotating neutron stars of q > 0.1. This implies 
that for rapidly rotating neutron stars, the effect due to 
the q n terms is significantly suppressed by that due to 



Q 2 . Thus, we conclude that the effect of the quadrupole 
is important in determining the ISCO even when q ~ 
0.1. On the other hand, the coefficients of q% and Q4 are 
smaller than those of q and Q^- Thus, their contributions 
are small. 

The formulas of E and I at the ISCO are shown in 
Appendix B. 

III. NUMERICAL STUDY 

To confirm the accuracy of the formulas derived in Sec. 
II, we compare them with numerical solutions. We nu- 
merically construct stationary, axisymmctric spacctimcs 
of relativistic rotating stars with a polytropic equation 
of state using Komatsu, Eriguchi, and Hachisu method 
Jl3| . Then, we estimate fhsco and -Risco as well as the 
multipolc moments of the numerically generated space- 
times. It is desirable to include all the multipole mo- 
ments up to Q4, and compare our analytic formulas with 
the numerical results. However, accurate numerical cal- 
culation of Q4 is difficult. In this section, we simply 
set Q4 = olQ\. For the Newtonian incompressible case 
(Maclaurin spheroid case), a — 15/7, and for Newtonian 
compressible case, we find a < 15/7. Due to general rel- 
ativistic effects, stars are more centrally condensed than 
those in the Newtonian case, so that we expect that a in 
general relativity is smaller than that in the Newtonian 
theory, i.e., < a < 2. Hence, we simply set a = 1. 
Fortunately, the coefficients of Q4 in Jlisco and flisco 
are small, so that our rough treatment does not affect 
the result much. In fact, we also set a = and 2, but 
our conclusion is not changed at all. 

A. Basic equations 

We consider the energy momentum tensor of an ideal 
fluid, 

= pb(l + s+ ^ju^u a +Pg IM ,, (3.1) 

where pt, P, e, it M , and g^ a are the baryon rest mass 
density, pressure, specific internal energy, four velocity, 
and spacetime metric. We adopt the polytropic equation 
of state, 

P = Kp\ = (r- l)p b e; r = l + i (3.2) 

n 
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where K and n are the polytropic constant and the poly- 
tropic index. 

We only consider the case when stars uniformly rotate 
around the z-axis. Hence, we set u r = u 6 = and u v — 
u°Q s , where fl s is the spin angular velocity of the rotating 
star. In this case, the fluid equations of motion are easily 
integrated to give |Q 

h 2 



(«°) 



2 = -h 2 (g u + 2g ttp ft s + g w tit) 



constant, 



(3.3) 



where h=l + e+ P/p b = 1 + Ts. 

Following Butterworth and Ipser fl5|l , we write the line 
element as 



ds 2 = 



3 2v dt 2 



r 2 sin 2 6B 2 e- 2u {dtp-u)dtf 



(dr 2 +r 2 d0 2 ), 



(3.4) 



where B, £, and to are field functions depending on r 
and 6. We present the field equations for these variables 
and briefly describe our numerical method in Appendix 
C. 

Once the field variables are computed, the ISCO on the 



equatorial plane is found making use of Eq. (2.8) together 



with Eqs. ( |2.5| )— (2.7). To examine the accuracy of our 
approximate analytic formulas derived in Sec. II, we also 
need to estimate the multipole moments. For M and q, 
we have the formulas as |15f| 



M 



d A x 



2h(u° 



'{l + e(2-r)> 



2r 2 sin^ 9(fL s - u>)ujB 6 (v^f e 



p b h 



(3.5) 



-^2 J d 3 xr 2 sin 2 6(fl s ~uj)B 3 (u a ) 2 e 



0n2„2C-4i/ 



Pbh. (3.6) 



For Q2 and 93, we estimate them by using the asymptotic 
behaviors of v and uj. For r — > 00, v, u> and B behaves 
as H 

M B Q M , Q 2 M 3 
r 



2qM 2 



3r 3 
6qM 3 



P2{y) + 0{r~% 
qM i 



(3.7) 



8 



35q 
M 2 



q 3 M q 



(5y'-l)+0(r- b ), 



B 



(3.8) 
(3.9) 



where y — cos£>, ^2(2/) = (3y 2 — l)/2, and B is a con- 
stant which is — M 2 /4 for spherical cases. When we nu- 
merically obtain v 1 to, and B, we can extract information 



of Q2 and qs at a large radius r ^> M near the outer 
numerical boundary (r = r max ) as 



<73 



M 3 
21r 5 



2 = 5 1TT / d V p 2i.y)v, 



16M 4 Jt] 



dy(5y 2 -l)(l-y 2 )uj. 



(3.10) 
(3.11) 



We evaluate Q2 and q% at various large radii, and find 
that Q2 quickly converges as r — > r max . This suggests 
that the error for estimation of Q2 is very small (we guess 
it is less than 1%). This is mainly because the coefficient 
of 0(r -4 ) part in v is not large. On the other hand, 
convergence of (73 at r — > r max is slow because the co- 
efficient of 0(r~ 6 ) part in to is large [^5|. To obtain q% 
accurately, it is desirable to attach the outer numerical 
boundary as r max 3> M, but to do that we need to take 
many grid points and hence need a long computational 
time. To save computational time, we use here an ex- 
trapolation method to estimate 93; i.e., we calculate (73 at 
several large radii which are not near the outer boundary 
(r ~ 3r max /4), and extrapolate true value of q^ at r — > 00 
by assuming that q^ behaves as qa(r) = 93(00) + C/r for 
the large radii, where C is a constant. Since this method 
is rough, we guess that the error of 93 may be ~ 10% in 
this method. However, the large error in 53 does not af- 
fect the following analysis much because the coefficients 
of q-j terms in fhsco and Aisco are fortunately small. 
The important quantities in our analysis are M, q and 
Q2. We mention that we have also evaluated M and q 
by using 



M 



dyv, 





3M 



8- 



3So 
M 2 



-1 -1 



3r 3 



AM 2 



x / dy{l-y 2 )u>, 



(3.12) 



(3.13) 



and confirmed that M calculated by Eqs. ( |3.5| ) and ( 3.1 2| ) 



and J calculated by Eqs. (3.6) and (3.13), respectively, 
agree very well. 

B. Results 

As the polytropic index, we set n = 1. In Fig. 1, 
we show the relation between p b $ {p b at r = 0) and M 
for the case n = 1. Note that in the figure, we plot 
non-dimensional quantities p b $K n and MK~ n / 2 . The 
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lower and upper solid lines denote the relations for the 
spherical star and for the rotating star at the mass shed 
limit, respectively. The filled squares denote data sets 
that are used for comparison with our analytic formulas. 
The dotted line is the critical line above which the ISCO 
ceases to exist. For sufficiently large p^o, stars are unsta- 
ble against radial gravitational collapse jll| . The dashed 
line divides the stable and unstable branches: The left- 
hand side of it is the stable region. (Here, for judging the 
stability against the gravitational collapse, we have ap- 
plied the turning point method shown in rcf. |p7[.) Since 
unstable stars are not realized in nature, we exclude nu- 
merical data in the unstable region. 

In Table I, we show all the data we use in our analysis. 
As pointed out in rcf. Q, Q 2 /q 2 is always greater than 
unity. (For the data in Table I, 2 < Q 2 /q 2 < 4. ) 

In Fig. 2, we show 

An (J) = 6V6M (^rsTo rical - ^isco(o) , 

Ar{l) = W (^sco liCal - ^co(o) , (3-14) 

as a function of Q 2 . Here, Oisco(0 and i?isco(0 (1 — 
I < 4) denote the analytic formulas in which we include 
terms up to 0(e l ). We also define A^ orr and Af cn ' in 
which we use the analytic relations for the Kerr met- 
ric. The open circles, crosses, filled circles, open triangles 
and open squares denote A^ orr , A a (l), A a (2), A a (3) and 
A a (4) (a = CI or r), respectively. 

From Fig. 2, we find apparently that the Kerr formulas 
for the ISCO are not appropriate at all. The formula 
-Risco(l) is as bad as the Kerr formula, but ^isco(l) is 
fairly good. The latter feature seems accidental: Since 
the relation Q 2 ~ aq 2 with a ~ 2 — 4 holds, cancellations 
between q 2 and Q 2 terms, between g 3 and qQ 2 terms, 
and among q 4 , q 2 Q 2 and Q4 terms occur. (In contrast, 
if Q 2 ~ q 2 , i?isco(l) w ih be bad, while Oigco(l) will be 
good. If Q 2 ~ 10q 2 , both formulas will not be good at 
all). 

It should be stressed that for small Q 2 , A^ crr and 
A^ crr (also A r (l)) are roughly proportional to Q 2 , while 
Aq (2) and A r (2) remain very small. Recalling that in 
the Kerr formulas f^igco and ^ISCO' the en?ec t of Q 2 as 
an independent variable is absent, this feature implies a 
substantial effect of Q 2 on the determination of the ISCO 
around neutron stars even for q ~ 0.1. 

We find the errors of the formulas f2isco(2) and 



-Risco(2) are always very small for small Q 2 < 0.2 
(q < 0.2 in this case). For large Q 2 > 0.2, however, 
the errors gradually increase. This indicates that the 
effect of the terms of 0(Q 2 ) (such as Q4) is not neg- 
ligible. Therefore, it should be necessary to correctly 
take into account the terms of 0(Q 2 ), in order to give 
an appropriate formula for very rapidly rotating neutron 
stars. Note, however, that Fig. 2 also indicates that 
A a (a = f2 or r) adequately converges to zero by adding 
higher order terms in e. Thus, unless q, Q 2 > 1, an ap- 
propriate formula for very rapidly rotating neutron stars 
of large q and Q 2 can be derived along the line presented 
in this paper. 

IV. SUMMARY 

In this paper, we have investigated the ISCO of a test 
particle moving on the equatorial plane around a rotat- 
ing object. We have derived fairly accurate approximate 
formulas for the angular velocity and circumferential ra- 
dius of the ISCO including mass, spin, quadrupolc, cur- 
rent octapole, mass 2 -pole moments of a rotating object. 
Our formulas show that the effect of quadrupole moment 
is important for determining the angular velocity and, 
in particular, the circumferential radius of the ISCO for 
rapidly rotating neutron stars of Q 2 > q |(|. 

In a recent paper, Miller, Lamb, and Cook B per- 
formed a numerical computation for analyzing the ISCO 
around rapidly rotating neutron stars, and pointed out 
that f^isco and in particular i?isco are not correctly de- 
termined by the first order analytic treatment in q. Our 
present study clarifies the reason for inaccuracy of the 
first order formula which is mainly due to the neglection 
of the quadrupole moment. 
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APPENDIX A: COEFFICIENTS FOR THE 
ERNST POTENTIAL 



1154583Q 4 2221047g 3 g\ 
739024 + 2956096 J ' 



(A2) 



FHP 

fl o,o 

a 0,5 
°0,6 



give 



M 20 M 



M 

163, a ,4 = Af 4 + 
.JAf 20 , M 30 M 
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21 

M 20 M 3 



+ 



5Af 20 M 2 4iM 30 J 8M 31 Af 



a 0,8 



a 0,9 



33 231 

6M 40 M 

+ ^l-< 

19iM 20 M 2 J 15M 30 M 3 



33 



33 



; 30 j 

143 

53Af 20 MJ 2 



429 

M 20 M 5 

143 + 3003 
36iM 3Q M 2 J 3M 3 



(M31 



a o,io 



143 13 
43iM 20 M 4 J 7M 30 M 5 

2431 
7M 20 M 7 



+ M 8 0(e 5 ), 

311Af 20 Af 2 M 2 
3003 

M40I + M 9 0(e 5 



221 

4423M 20 M 4 M 



M 10 O(e 5 ), 

202M 20 M 3 J 2 



4199 
462iM 30 M 4 J 
4~L99 
+ M n O(e 5 ), 



138567 
28M" 



12597 



+ ^2T( M31+M40 ) 



(Al) 



where M 20 = AfAf 2 + J 2 , M 30 = i(S 3 M - M 2 J), M 31 = 
-5 3 J - Mi, and Af 40 = M 4 M + S 3 J §. 

Using Eq. flgjj ), 0^,0 (1 < 3 < 5) and a 2i ,i (1 < j < 
4) up to 0(e 4 ) are calculated as 



a 2 ,o 



«4,0 



06,0 



M 3 

-— (l-Q 2 ), 



M 5 
"56" 



(21 + 10 9 2 -52Q 2 + 2ig 4 ), 
313g 2 4717Q 2 



M — — 
16 

?2 



924 



289q 2 Q 2 38Qj 



3696 
15ig 4 , 17g 3 g 



^8,0 



M 



1848 
35 

128 + 
32171<? 4 
384384 
5571Q4 



77 176 66 
31303g 2 101373Q 2 
64064 64064 
22513g 2 Q 2 57535Q| 



27456 
3343g 3 g\ 



aio.o = M 



4576 6864 J ' 
63 12974415g 2 
256 ~ 20692672 
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34944 



153744405Q 2 
82770688 



6074105g 4 1686890165g 2 Q 2 2883979925Q 2 , 



22284416 



869092224 



869092224 



and 



a 2 .i 



3iM 4 



-9 + 93), 



a 4 ,i = iM I — 



15g , 15g 3 29gQ 2 13q 3 



28 28 A P 



0/ 35g 3763g 3 95213gQ 2 30743g 3 

a« 1 = iAz 1 1 . . 

V 16 2002 24024 6864 J' 

,, 10 /3l5g 579957<j 
08,1 = iM f 



128 155584 

154403<?Q 2 876263<j 3 
19448 155584 



(A3) 



APPENDIX B: FORMULAS OF E AND £ AT THE 

ISCO 



From Eqs. (2J3) and (2/7) with an approximate for- 
mula of /Cisco shown in Sec. II, E and I at the ISCO are 
derived as 

£isco = 0.94281 - 0.03208g - 0.02975g 2 + 0.00794Q 2 

- 0.0341g 3 + 0.0198gQ 2 - 0.0019g 3 

- 0.0440g 4 + 0.0404g 2 Q 2 - 0.0033Q^ 

+ 0.0005Q4 - 0.0062gg 3 , (Bl) 



Cisco 
M 



3.4641 - 0.9428g - O.Uteq 2 + 0.1879Q 2 

- 0.3953g 3 + 0.2996gQ 2 - 0.0392g 3 

- 0.4470g 4 + 0.4944g 2 Q 2 - 0.0505Q^ 

+ 0.0093Q 4 - 0.0926gg 3 . (B2) 



On the other hand, exact solutions for the Kerr case (U 
are expanded as 



Ef^ = 0.94281 - 0.03208g- 0.02182g 2 
- 0.01633g 3 - 0.01294g 4 

pKcrv 

isco = 3.46410 - 0.94281g- 0.25660g 2 



M 



0.13531g 3 - 0.08791 9 4 . 



(B3) 



(B4) 



Hence, as in the case of f^isco an d -RiscO; the error of 
coefficients is less than 10 -4 for 0(e 2 ) terms, ~ 10~ 3 for 
0(e 3 ) terms, and ~ 10~ 2 for 0(e 4 ) terms. 

If accretion of matter occurs from the ISCO to the 
central rotating body, the value of q of the central star 
increases as 



8 



1 M\ M 



ZqEisco , 



(B5) 



where p is the rest mass of the accreting matter. Using 



Eqs. (Bl) and (B2), an approximate formula for 5q is 



given as 



Sq= M ( v 3 - 4641 ~ 2 - 8284 ? ~ 0.3802g 2 + 0.1879Q 2 

- 0.3358g 3 + 0.2837gQ 2 - 0.0392g 3 

- 0.3788<? 4 + 0.4548g 2 Q 2 - 0.0505Q^ 



0.0093Q 4 - 0.0887^3 



(B6) 



Thus we obtain the following conclusions from Eqs. ( Bl ) , 
If accretion of matter occurs from the 



(|B2D, and 

ISCO, the effect of Q2 are that (1) the maximum energy 
release efficiency of accreting matter slightly decreases, 
and that (2) the rate of increase in the angular momen- 
tum and Sq of central body is increased. 



APPENDIX C: FIELD EQUATIONS FOR 
COMPUTING ROTATING STARS AND THE 
NUMERICAL METHOD 

In this appendix, we show field equation for v, u>, B, 
and £ for obtaining axially symmetric rotating stars, and 
briefly mentioned our numerical method to solve them. 
The field equations are as follows |l5| . 

2 1 

V,rr + ~V,r + ^{"-.vyi 1 ™ J/ 2 ) ~ ^V v ,y\ 

= 4ire 2<: p b [2h(u ) 2 - {1 + (2 - T)e}e~ 2 "} 
+ 1 -e-^B 2 (l-y 2 ){rW r + ul(l-y 2 )} 



(CI) 



= -I67re 2C p 6 /i(u°) 2 (l -u>)+ Uv, r - 

1 - y 2 ( Av ^ 



35,, 



r 2 V " " r> 



(C2) 



B,rr + -B, r + \{B, yy (l - y 2 ) - 3yBj 

= l6irBe 2C - 2,/ P, (C3) 
l v = [{B + B, r r) 2 (l - y 2 ) + {By - B, y (l - y 2 )} 2 }- 1 

[(B + B, r r){-e- 4l/ Lu, y Lu, r r 3 B 3 (l - y 2 ) 2 /2 
+ rB >r y(l - y 2 ) - B, r ry + 2v v v, r rB(l - y 2 )} 

\{By - 5,(1 - y 2 )}{e-^B\ 2 (l - y 2 )(u 2 r r 2 



w 2 v (l - y 2 ))/2 - 5 rr r 2 - B <r r - 2Bv 2 r r 2 



,y 

B )1W (1 y 2 ) 3yB, y 



2^ 2 (l-y 2 )S}], (C4) 



where instead of u>, we introduce uj = ujflj 1 because uj oc 
£l s for slowly rotating cases. 

We solve Poisson type equations for uj (uj), and B 
as the boundary value problem. That is, we solve these 
equations imposing boundary conditions (3.7) — (|3~^) at 
r — r max > 10M , regularity condition at r = and 
6 = 0, and reflection symmetry conditions at 6 = ir/2. 
On the other hand, the equation for £ is solved using sec- 
ond order Runge-Kutta method imposing the boundary 
condition at0 = Oas£ = lri-E? . 

Using field variables, Eq. (|3.3[) is rewritten as 



(1 + nKp\ ln ) 2 [-e 2v + tt 2 (l - uj) 2 B 2 e 
— C(= constant). 



-2v 1 ■ 2 , 

r sm 1 



(C5) 



Here, C , K, and Q s are constants determined in iteration 
which is carried out in the following manner |l3| : 

1. First of all, we fix the coordinate radii of stellar 
surfaces at equator and pole. Also, we fix the value 
of pb at r = 0. 

2. We give a trial density configuration for p b . 



3. We solve field equations (C1)-(C4). 



4. Since we have constraints imposed at procedure 1, 
we can determine C , K, and il s using Eq. (|C5|). 



5. Using Eq. (CS) with C, K, and Q s determined at 
procedure 4, we calculate a new trial configuration 
for p b . 

6. Return to 2, and repeat procedures 2—5 until a 
sufficient convergence is achived. 

Numerical computation is typically performed taking 
(N r , N e ) = (800, 80) grid points, where N r and N e denote 
grid points for < r < r max and < 6 < n/2. We use 
homogeneous grids in r and y = cos 6. We check our 
numerical code by comparing our numerical results for 
non-dimensional quantities MK~ n l 2 , JK~ n , RK~ n l 2 , 
and Q s K n / 2 with another results presented previously in 
ref. jl6| for the case n = 1. We find that our numerical 
results agree well with those in ref. |l(| (inconsistency is 
less than 1%). 
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Table I. Numerical data for n = 1 polytrope which 
we adopt in Sec. III. Note that R e denotes the circum- 
ferential radius at the stellar surface on the equatorial 
plane. Because MO s (e \) is an invariant quantity irre- 
spective of scaling of if, we can calculate the spin angular 
frequency of stars from the following data approximately 
as 37OHz(x/O.O16)(1.4M /M). 



Pb.oK" 


MR- 11 ' 2 


R e K- n l 2 


n s K n / 2 


q 


Q 2 


93 


i?isco/6M 


ttisco6V6M 


.200 


.158 


.870 


.061 


9.48E-02 


3.53E-02 


6.9E-03 


.953 


1.068 


.217 


.160 


.852 


.064 


9.43E-02 


3.19E-02 


6.1E-03 


.953 


1.069 


.236 


.162 


.833 


.067 


9.38E-02 


2.91E-02 


5.4E-03 


.953 


1.069 


.257 


.163 


.814 


.070 


9.34E-02 


2.66E-02 


4.8E-03 


.952 


1.070 


.280 


.164 


.795 


.073 


9.31E-02 


2.41E-02 


4.4E-03 


.952 


1.070 


.305 


.164 


.776 


.076 


9.28E-02 


2.22E-02 


4.0E-03 


.952 


1.070 


.206 


.160 


.869 


.098 


1.50E-01 


8.45E-02 


2.6E-02 


.931 


1.105 


.224 


.162 


.851 


.103 


1.49E-01 


7.66E-02 


2.3E-02 


.930 


1.106 


.243 


.163 


.832 


.107 


1.49E-01 


6.99E-02 


2.1E-02 


.929 


1.108 


.265 


.164 


.813 


.112 


1.48E-01 


6.36E-02 


1.8E-02 


.928 


1.109 


.289 


.165 


.794 


.116 


1.47E-01 


5.82E-02 


1.7E-02 


.927 


1.110 


.314 


.165 


.775 


.121 


1.47E-01 


5.36E-02 


1.5E-02 


.927 


1.111 


.222 


.163 


.863 


.144 


2.12E-01 


1.53E-01 


6.6E-02 


.908 


1.145 


.241 


.165 


.844 


.150 


2.11E-01 


1.39E-01 


5.9E-02 


.906 


1.149 


.262 


.166 


.825 


.156 


2.10E-01 


1.27E-01 


5.3E-02 


.904 


1.152 


.285 


.166 


.806 


.163 


2.10E-01 


1.17E-01 


4.8E-02 


.902 


1.154 


.310 


.167 


.787 


.170 


2.09E-01 


1.07E-01 


4.3E-02 


.901 


1.157 


.237 


.166 


.859 


.181 


2.60E-01 


2.11E-01 


1.1E-01 


.891 


1.176 


.257 


.167 


.840 


.189 


2.59E-01 


1.93E-01 


9.9E-02 


.888 


1.181 


.280 


.168 


.820 


.197 


2.58E-01 


1.77E-01 


8.9E-02 


.885 


1.186 


.304 


.168 


.801 


.205 


2.57E-01 


1.63E-01 


8.1E-02 


.883 


1.190 


.231 


.167 


.876 


.206 


3.03E-01 


2.86E-01 


1.7E-01 


.882 


1.196 


.251 


.168 


.857 


.214 


3.01E-01 


2.61E-01 


1.6E-01 


.877 


1.203 


.272 


.169 


.838 


.223 


3.00E-01 


2.40E-01 


1.4E-01 


.873 


1.210 


.295 


.170 


.818 


.232 


2.99E-01 


2.21E-01 


1.3E-01 


.870 


1.216 


.249 


.172 


.883 


.258 


3.71E-01 


3.84E-01 


2.9E-01 


.862 


1.237 


.270 


.173 


.863 


.269 


3.70E-01 


3.54E-01 


2.6E-01 


.856 


1.247 


.292 


.173 


.844 


.279 


3.68E-01 


3.27E-01 


2.4E-01 


.851 


1.257 


.274 


.176 


.887 


.307 


4.28E-01 


4.53E-01 


3.9E-01 


.844 


1.276 


.294 


.177 


.869 


.318 


4.26E-01 


4.22E-01 


3.5E-01 


.837 


1.288 
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FIG. 1. Relation between mass (MK' n/2 ) and density at 
r — (pb,oK n ) for n=l polytrope. Lower and upper solid 
lines denote the relations for the spherical star and for the 
rotating star at the mass shed limit. Dotted line (with open 
circles) divides two regions where the ISCO exists or not: 
It exists below the line. Dashed line divides the stable and 
unstable branches: The left-hand side is the stable region. 
Filled squares denote data sets which are compared to our 
analytic formulas. 



FIG. 2. An and A r as functions of Q2 for n = 1 poly- 
trope. Open circles, crosses, filled circles, open triangles 
and open squares denote A* orr , A«(l), A„(2), A„(3) and 
A a (4) (a — Q or r), respectively. 
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